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ABSTRACT 

We present the results of a numerical, 3-dimensionai study 
of cc'.vection in a self-gravitating sphere of Boussinesq fluid 
with a Rayleigh number of 10^*^ and a Prandtl of 1. The velocity 
and temperature are computed by using spectral methods (spherical 
harmonics, in the horizontal and finite -differencing 

in the radial directions. An eddy viscosity and diffusivity 
are needed to model the sub-resolution flow. The amplitudes 
of the eddy viscosity and diffusivity are determined by requiring 
that the time-averaged kinetic and thermal energy spectra remain 
unchanged when the limit of spatial resolution varies . For Ray- 
leigh numbers much less than 10^*^ the flows do not have well- 
defined inertial ranges and an eddy viscosity and diffusivity 
cannot be assigned in this self-consistent manner. By computing 
the energy spectra as well as the detailed energy budgets as a function 
of wavenumber, we show that for Rs = 10^*^ there is an inertial range 
for the modes corresponding to spherical harmonics with £>6 . We 
also show that the velocity field becomes nearly isotropic at small 
wavelengths. Intermittent bursts of convective flux and energy prop- 
agate inward from the outer boundary-layer. The bursts are shown 
to cascade from large to small wavelengths . 

Subject headings : convection - stars: interiors - 

hydrodynamics 




I. INTRODUCTION 


This is the third in a series of papers on convection in 
spheres. In it we extend our numerical calculations to high Ray- 
leigh numbers. In the secoxid paper of this series, the Rey- 
nolds and Peclet numbers were small enough so that the small length- 
scales where dissipation occurs could be resolved. In this paper 
the dissipative lengthscales are much smaller than the numerical 
reso].ution, so it is necessary to model the small wavelength end 
of the kinetic and thermal energy spectra. Clearly, no attempt 
to compute the convective flows in the cores of stars (where the 
Reynolds numbers are greater than 10^*^) is able to resolve 
all of the scales of motion so it is necessary to model the 
unresolved part of the flow. For example, the mixing-length treat- 
ment of convection and the spectral method proposed by Ledoux et al . 
(1961) both assume that buoyancy drives the convection in a narrow 
band of large wavelength modes. Both methods utilize the 

model of homogeneous, isotropic turbulence in which buoyancy is 
not important to compute the effects of the small-scale unresolved 
motions . 

In an actual sphere of convecting fluid we might expect that 
the buoyancy effects are felt over a wide range of wavelengths, and 
that the fluid which is strongly driven by the buoyancy in the 
radial direction is neither uniform nor isotropic for an even 
wider range of wavelengths . In this paper we explicitly compute 
the fluid's velocity and temperature at the largest wavelengths 
and confirm that the fluid is anisotropic and that 
the energy spectra are not smooth functions of wavelength but 


have a large amount of fine structure. At smaller wavelengths 

we also compute the flow and show that it becomes 

isotropic and that the energy spectra develop inertial ranges 
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with the smooth Jc power law associated with homogeneous 
turbulence. We are then able to model the effects of the small 
scale fluid motions by assuming that the calculated intertial range 
smoothly continues beyond our limit of spatial resolution. In § 2 
we show how to parameterize the transport properties of the unre- 
solvable inertial subrange with eddy viscosities and diffusivities . 
The time-averaged results of our calculations are presented in §3. 
In the fourth section we examine the time-dependent fluctuations 
in the energy spectra and show how they cascade fro-m large to small 
wavelengths. A summary of our conclusions appears in §5. 
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II. EDDY VISCOSITY AND DIFFUSIVITY 
A) Motivation 

Unfortunately, it is impossible to compute numerical solu- 
tions to the Navier-Stckes equation for large Reynolds number flows. 

12 

The Reynolds number in the convective core of an 0 star is 10 , 

which means that the ratio of the pressure scale-height to the 

g 

dissipative length-scale is 'v 10 . Unless a numerical method 
has the capability of resolving this range of scales, the resulting 
solution will lack some essential physics. To circumvent their inabil- 
ity to resolve the dissipative scales , hydrodynamicists ' typically 
model the dissipation by replacing the true kinetic viscosity, 

V, with an eddy viscosity, v , with v^»v. It is usually argued (c.f. 
Deardorff, 1972 and references therein) from dimensional analysis 

that V should be of the form 
e 


V = c L V 
e res res 


( 2 . 1 ) 


where L and V are the length and velocity at the limit of 
the spatial resolution and c is a constant of order unity that 
is independent of both position and . Similar reasoning shows 

that an eddy diffusivity, of the same form as equation (2.1) 

should be included in the heat diffusion equation. More compli- 
cated models for eddy viscosities exist (c.f. Kraichnan, 1976). 
However, a simple test of the validity of an eddy viscosity of the 
form of equation (2.1) is to see if we can find a value of c (inde- 
pendent of a position) such that all solutions computed with this 
value of c, but with different limits of spatial resolution, have 
the same energy spectra. This type of test is easily performed 










with a modal analysis. 


B) Definitions and Their Implications 


Let the temperature T be a sum of its resolvable T and unresolv- 
able T' parts. If we assume that the radial resolution is much 


greater than the angular resolution, then we may write 
^cutoff ^cutoff 

^ t. T. = I T, 

il = 0 .. ^ = 0 


( 2 . 2 ) 


m,Y 


and 




£ = ~^. + l 

cutof r 




z 


m,Y 


Z = Z , -^+1 

cutoff ( 2 . 3 ) 


The sums in equation (2.2) - (2.3) are over with y equal to 

R and I, where 


Y^’^’m = 2(2^)^^^^ Re (Y^’”^) m 0 


2(tt)^^^ Y^’° 




m 


= 0 


(2.4) 


Yl,5.,m ^ 


2(27r)^'^^ Ira m ?? 0 


(2.5) 


m = 0 


,5. ,m, 


and where Re (Y^) and Im(Y'^’'“) are the real and imaginary parts of 
the spherical harmonic. We use the notation throughout that any 
quantity written with subscripts y,Z,m is a function only of 
radius and time and corresponds to the (Yj^^jIh) component of a 
quantity 5 i . e . , 




f Y~<’ 

4tt 


£ ,m 


f (r , 6 ,(|) ,t) dS^ 


( 2 . 6 ) 


Writing v=v + v’, we can define an eddy diffusivity each 

^^cutoff 


V • 7C? „ ^ - [(v-V) T]^ p yY,^-,m 

Y 5 X. ,m — Y , £ ,m 


= - [(v-V) T . ^ 

- Yj^»m 


(2.7) 


All of the quantities on the left-hand side of equation (2.7) are 

numerically calculable except for If v and T were known 

over all scales, we could then solve for 'R ( 5 , ^-,r,t) and 

le cutoff’ ’ 

determine how strongly it depends on Y)^?m»t,t, and 

Similarly we can define two eddy viscosities for each (YjJijm) compon- 
ent of the velocity. There are two viscosities because the 
velocity has two degrees of freedom. 

The toroidal eddy viscosity, ,t) , is defined: 

,{7x[V. (v„^’^’^Vv3. 

( 2 . 8 ) 


e • {V X [(v*V)v - (v-V)v]} 
r — — 


=. e 


Y,^,m ~ — ''i("Y,A,m 


and the poloidal eddy viscosity, v_^ ’ ^ ^cutoff 


{C(vV)v - (vV)v] - 0 = e • [V-(v ^'’-^’^Vv)]) „ 


(2.9) 

The pressure, that appears in equation (2.9) is the 

part of the resolvable pressure that is due to nonlinear terms 
involving v' : 


P m Y^’^’^)= -{V- [(vV)v - (J-V)J]}^ j, ^ 


( 2 . 10 ) 

It is necessary to include P^^ in equation (2.8) to make 
[(v»V)v - (v*V)v - VP^^] divergenceless. The resolvable pressure 

/S 

P p is then equal to 
y 5 J6 5 In 






sim 


P = p'-*- + p- 

YjAjm 


( 2 . 11 ) 


where 


^2(pex ^ -{V- C (J-V) v) n 

Y,^,m - - y,H,m 


+ r"^R_P„ 8(r^T., „ _)/ar (2.12) 


s r 




Equations (2.8) and (2.9) provide a complete description for the 

A • 

resolvable part of C(v*v)v - (v*V)v - Vp 


/\ A 


[(v.y)v - (v.v)v - VP^"^] = 


^ ^ CA(i-+l)]"^ {VaCr^e -V(v Y,^',!!!^^) ]/ar 

Y>^jm r p - 


- r 


V^Cr^’e *V(v _ 

, r p — r 


/S /\ 

e 


X V Cr^e^- {Vx [ V ( ^ ’^Vy) ] > 


(2.13) 


T p 


If V =v one might conclude from equations ('’.8) and (2.9') that 


[(v*V)v - (v*V)v - VP^^] is equal to V'(VpVvJ, but this is not 


p _■ 


true because the latter quantity is not divergenceless. When modeling 


v' with an eddy viscosity, one must be careful that the model 

2 


is realizable; that is (v' ) should never be negative. This problem 
is frequently encountered (c.f. Schlesinger, 1978) but does not 
occur in our modal definitions of eddy viscosity. 

For the eddy viscosities and diffusivity to' be numerically 
useful tools, they must be independent of y and m. At large wave- 
numbers, we expect the temperature and velocity to be isotropic. 
Furthermore, if the nonlinear terms are dominated by local inter- 


actions in phase space {i.e. the interaction of v , with 

Y ) ^ 5^1 


T „ „ „ is an important contribution to C(v*V)T] „ only if 

YsYs^l — 










then we expect 4 ^ to be independent of y and m for 


large values of 1. For small values of I the flow is not isotropic 


0 Y 0 m 

and ’ ’ will not be independent of y,m. However, if -j-of-^ suffi- 


ciently large, then 


^ . „I«KCv.V)T]^ „ for il<<£ (2.14) 

'e Y,x,m‘ ' — Y,x,m' cutoff 


Therefore , small errors in eddy diffusivity caused by assuming that 
is independent of y and m for low values of H will be negligible in 
comparison to the explicitly calculated nonlinear advective terms. 


Assuming that the eddy diffusivity and viscosities are indepen- 
dent of Y and m, we parameterize them as : 


J = 2tt[ 


(A 


cutoff cutoff 


-+1)]”^^^[KE(£ 


cutoff 




K 

(2.15) 


^p ■ ^^^^cutoff ^^cutoff'*'^^^ ' ^-^^^^cutoff ^p 


(2.16) 


"t = 2^'^^cutoff^^cutoff"l^''^'''f^^E^Woff ’^^/27rr2]l/2 


(2.17) 


where KE(£,r) is the kinetic energy of all modes with wavenumber I at radius r (see 
Marcus 1980a, hereafter referred to as Paper II). The length and velocity at the 


resolution limit are taken to be ^irC £^^^^^^+1) ] and 


,r) /2iTr^ respectively. The eddy coefficients 


C Cm, and C, are functions of I, as well as r and t. However, 

P 5 i K 

if lies in the self-similar inertial range, dimensional 

arguments can be used to show that Cj^, and 

should be independent of i . The gist of the argument is that 
the rate at which energy is transferred from modes of wavenumber 
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% to modes with wavenumber greater than il , is (Tennekees 

cutori 

and Lumley 1972) 


SCJO" ^.)/S(£ ..) 

cutoxt cutorr 


(2.18) 


where S(^-) is the strain and E(il) is the energy associated with 

modes of wavenumber ii.If then and so the 

. ' . 4/3 

rate is proportional to ^ eddy viscosity in the 

modal Navier-Stokes equation, this same dissipation rate is proportional to 


E(£) V = C ( 

P P 


£/£ ^ -.) 
cutorf 


4/3 


(2.19) 


if Iv Setting the two expressions for the rates equal, we see that Cp 
is independent of %. 


C) Self-Consistency and a Method of Solution 


The most difficult task in implementing the eddy viscosity 

is choosing the values of C^, and C^. They are usually chosen 

(Deardorff, 1972; Siggia, 1978) by varying their values until the 

-5/3 

resulting kinetic and thermal energy spectra have smooth, k ' 
power laws. This method is not only arbitrary but somewhat expen- 
sive in that several solutions with differing values of the eddy 
coefficients must be computed before the correct values of the 
eddy coefficients can be determined. 

We shall establish the values of the eddy coefficients by 
appealing to the self-consistency requirement that , Crp, and 
Cv are indecendent of £ , — Snecif ically, the time-aVeraged 

thermal energy , poloidal kinetic energy, and toroidal kinetic 

enerc-y stectra should not change when £ , is changed 

oj ^ = cutorr 




10 - 


•from 12 to 10. A necessary condition for the thermal 

energy spectra to be the same is that the rate at which the ther.- 
mal energy cascades in and out of the modes with wavenumber 
i due to tho advective term is independent of The rate 


of this energy cascade is ^ 


Y,m 


Tv 0 (, . For the cascade rate 

Y5JO5JU — y^x^jin 


with to be equal to the rate with requires that 

e e y ,1 ,m 




+ [v +v 1 *v Ft YT 

^-11 -12^ ^ni ■^12^'^Y,^.m 


( 2 . 20 ) 


where T 5 ^ jniQ2) and /? T5^>niQg) calculated from eq.(2.15) 

A A 

using ^au+off to 12 and 10 respectively and where T and v 

are the resolvable temperature and velocity ^g^•toff "TO . 

Using equation (2.15) for-rl in equation (2.20) we obtain 


a linear differential equation for Cj^ (r) which may be numerically 
solved: 


V 


'K 




CA(r) ^ T t). 

Y,ni 

Y,m 


( 2 . 21 ) 


Y,m 


where 


and 


A(r) S (2TT)^^^r"^[(12*13)"^'^^KE(12,r)^^^-(10-ll)KE(10,r)^'^^] 

( 2 . 22 ) 

where . is the 2-dimensional Laplacian, X^p= (V^-l/r9^/8r^r) . 
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Now, obtained from equation (2.21) will obviously be a func- 
tion of r. Also, we can obtain different values of by choosing 
a different value of Z in equation (2.21) or by choosing a different 
pair of values for is a test of the self-consistency 

of our eddy diffusivity model if we can show that does not 
vary much with radius or depend strongly on our choices of Z and 
^cutoff* examining the kinetic energy budgets of the poloidal 
and toroidal velocity, we obtain differential equations similar 
to (2.21) for C and Crn* 

An advantage of using equation (2.21) to compute at each 

timestep is that the iterative values converse rapidly* To see the raoid 
convergence, let us assume that we have underestimated C^. A 
small value of makes the thermal energy extraordinarily large 
at; the il=ll and Jl=12 wavenumbers since it is not dissipated and 
piles up near ^^^.^off Marcus 1980b). The upward curl of the 

energy spectr'um makes the right-hand side of equation (2.21) large 
and thereby increases C^. The convergence of the equations that 
'govern and C,p is even more rapid because it can be shown that 
Cp and are nearly proportional to 1/A(r) or 
l/{CKE(10,r)]^'^^-. 8397CKE(12.r)]^''^} . If or is underesti- 
mated, the upward curl on the energy spectrum at i=12 causes this 

denominator to shrink, causing Cp and C,p to increase rapidly. L5kewiss, 
if Cp is too large, KE( 12 ,r) dissipates to zero, the denominator 

becomes large, and decreases. 

D) Numerical Results 

We have integrated the modal equations of convection (Paper II) 
with the eddy viscosities and conductivity described in this section. 
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0 

For a Prandtl numbei- of 1 and a Rayleigh number less than 10 , 
the values obtained for the eddy coefficients C^Cr) , C^(r) and 
Cj^(r) varied widely with radius, Z and pair of ^cutoff's' 
have concluded that our model of the eddy viscosity is not valid 
in this regime. This result should not be surprising since the 
Reynolds numbers are less than 500 and the high wavenumber modes 
are neither isotropic nor lie in an inertial range. 

For a Rayleigh number of 10^^ the eddy coefficients do con- 
verge. With "to 10 and 12 and Z equal to 8 in 

equation (2.20), the value of changes by less than 15% from 
r=0 to r=0.8. The time- averaged value of averaged over 
o^r^. 8 is 0.52. In the turbulent boundary-layer between r=0.8 
and r=1.0,Cj^ is 20% lower than this value. The coefficients 
and Cp are similarly well behaved with average values of 0,55 
and 0.46 respectively. Using the velocity and temperature pro- 
duced from these eddy coefficients in eq. (2.21) with £=6 and 
^cutoff 8, we have computed a new averaged value 

of that is ^ 10% lower than the old value. The de- 

crease in eddy diffusivity is due to the fact that the calculated 
2. = 6 energies are higher (due to production) than we would have 

predicted by using the value of the thermal energy at £ = 12 
■ -5/3 

and extrapolating to £ = 6 using a k spectrum (see §3). 

New average values of and were also computed based on 
£ = 6 and = 10 0* Tbe value of Cp was 10% lower, 

but was 10% higher than its original value because the 
£ - 6 toroidal kinetic energy is lower than we would have 
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predicxed (assuming a k 


-S/3 


spectrum and knowing t>Ve value of 


the toroidal kinetic energy at A = 12). The toroidal kinetic 
energy is relatively low at low wavenumbers because it is not 
produced at the low wavenumbers by buoyancy, it is produced 
at the high wavenumbers by isotropisation of the velocity. 

The new values of Cp, C,p , and like the old values, 

are only weakly C'vl2%) dependent upon radius. We have computed 
the time-averaged eddy coefficients with 3 other possible combi- 
nations of IL and and find only a weak (''>15%) dependence 

on these values. Because the eddy coefficients are nearly indepen- 


dent 


of I, r, and the pair of ^c^.toff ' s’ conclude that our choices 


of the functional forms of the eddy viscosities and diffusivity 

are self-consistent for R =10^°. It is also reassuring that the three 

s 

eddy coefficients are nearly the same and -are all of order unity. Tor 

the remainder of the numerical results presented in> this paper , 

we have held C^, C,p and fixed at their averaged values of 0.46, 

* 

0.55, and 0.52 respectively. 
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III. TIME-AVERAGED NUMERICAL RESULTS 
A) Temperature and Flux 


We have computed a numerical solution for convection in a 


ao 


sphere with a Rayleigh number, R , of 10 , and a Prandtl number, P , 

s ^ 


of 1, and with the boundary conditions and dimensionless units 
described in Paper II. , The heat source, H(r), is constant for 
r<0.3 and zero outside of this radius. We have used the time- 
averaged eddy coefficients described in §2. Although the con- 
vective solution is time-dependent, the energy spectra are time- 
independent when averaged over several of the largest eddy turn- 
over times; the results presented in this section have been 
averaged over 'v 30 eddy turnover times. We use the notation that 
single brackets, < >, around a quantity indicate an average over 
d,(p and double brackets, << >>, indicate an average over B,(p, and t. 
The horizontal and time-averaged temperature as a function of 
radius, <<T(r)>>, is shown in Figure 1. The fluid is nearly iso- 


thermal except at the outer boundary which has a thickness of 


g T — X / 3 

'V 0.03. By an asymptotic expansion it can be shown that 1^1 

_4 

=5x10 faraway from the boundary layer. The time-averaged numer- 


ical results are not inconsistent with this value, but at any par- 

3T 


ticular value of radius and time the fluctuation in |-^| can be 


10-100 times this value (see §3). The central temperature is 


<<T( 0) >> = 0 . 036, which should be compared to the central temperature 


■of 0.68 calculated with R =10^, P = 10 (Paper il) and a central 

s r 


temperature of 4.0 for a fluid in conductive equilibrium. Although 

3T 


3 T —1 / 3 

•^1 scales as j the central temperature scales as |-^1 x 


boundary-layer thdckness. Since the turbulent boundary-layer thickness 
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has a slight Rayleigh nturiber dependence , <<T(0)>> scales only 

■"1/ 3 

approximately as (R ) . Figure 1 shows a slight bumpatr=0.3 

* s 

which is the boundary of the heat source. This bump is not nearly 

4 

as Drominent as it was for R =10 . and we conclude that at large 
* s ' ° 

Rayleigh numbers, the convective velocity is very efficient in 
smoothing out the temperature perturbations due to the inhomogeneities in the 
heat source. The mean temperature curve also shows that the 
temperature gradient, << ^ is positive near the origin and 

near the outer boundary layer. In these two regions, the conduc- 
tive flux transports energy toward the center of the star; the 

time-averaged convective flux due to all modes, Total F , is 

con 

greater than the steady-state flux of the heat source, 

F^_(r) = (47rr^)~^ /^H(r')dr'. We have plotted the ratio of these two 
s s o 

fluxes in Figure 2. Far away from the boundary-layer this ratio 

-1/ 3 

should be l-(R ) . Convective overshoot is indicated where 

s 

Total F is greater than 1. Because Total F ~F ^ and is 

con ss ° con ss 

so much greater than the conductive flux, there can be large (10- 
100) fluctuations in the conductive flux (i.e. temperature grad- 
ient) without affecting the total flux. Plotted on the curve in 
Figure 2 with error bars is the variance of Total F /F^„ (The 

variance a of a quantity x is with respect to time, 

■ 2 1 / 

a = [<<(x-<<x>>) >>] “). The variance of the total (convective + 

conductive) flux divided by F is nearly identical to the variance 
shown in Figure 2. The increase in a with radius is due to the 
fact that the local Rayleigh number in the surface boundary- 


layer is near its critical value. Secondary convective insta- 


bilities in the boundary~layer cause time-dependent bursts in the 
convective flux. These fluctuations propagate inward and diminish 
in amplitude as they approach the origin. Notice that the variance 
at the surface is 18% of the total flux. 

B) 2- and 3-Dimensional Spectra 
The time-averaged kinetic and fluctuating thermal energy spec- 
tra (see Paper II for definitions) at r=0.5 as functions of the 
2-dimensional wavenumber, I, are shown in Figure 3. Both spec- 
tra behave as power laws at large 1, but for £-6 the kinetic energy 
spectrum is nearly exponential and obviously shows the effects 
of production. The thermal energy spectrum has a noticeable "cliff" 
between and Z=S. Each spectrum is normalised by the total energy 
at r=0.5 (due to all modes). 

• . 2/3 

In dxmensionless units the kinetic energy (which scales as ) 

at r = 0.5 is ^ KE ( Jl ,r=0 . 5 ) = 6 . 26xlo'^ . Approximately 91% of tliis 

energy is in the £=1 modes. The fluctuation thermal energy (which 

scales as R ^^^) at r=0.5 is '4.52x10”® and approximately 35% of 
s 

this energy is in the £ - 1 modes. The variances of the H = 1 
component of the kinetic and thermal energies are .08 KE(1, 0.5) 
and 0.2 TE(1, 0.5), respectively. The energies of the three 
i - 1 modes and completely change in one 

eddy turnover time. If the three modes were statistically inde- 
pendent , then the X- = 1 kinetic and thermal energy would each 
have a fractional variance of 1//3” . Because our numerically 

determined variances are much smaller, the three 2, = 1 components 
of the kinetic and thermal energies must be very well correlated 
with each other. The 2-dimensional kinetic energy spectrum has been 
deconvolved into a 3-dimensional spectrum and plotted as a function 
of the log of the 3-dimensional wavenumber, k, in Figure 4. The 
dashed line is k and is provided for comparison with a Kolmogorov 
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spectrum. The spectrum clearly shows a production peak at the low 
wavenumbers, a cliff at k^3 and an equilibrium range for k>6. The 
wavenumber of the Kolmogorov length calculated from rate of kinetic 
energy dissipation (§3.d) is '\^3xl0 . As discusvsed in Paper II, approxi- 
mations are introduced into the transformation of a 2-dimensional 
spectrum into a 3-dimensional spectrum, so. the 5/3-power-law behavior 

of Figure 4 should be viewed with some caution. 

The kinetic energy KE(£,r) as a function of radius 
is graphed for each value of S, in Figure 5. The curves are qual- 
itatively the same as they were for = j (Paper II). 

Except for Jl = l and Ji = 12 each curve is aero at the origin, rises to 
a peak at 0.5<r<0.8, decreases to a local minimum, and then increase 
sharply in the outer boundary-layer. The peak shifts to larger 
values of radius with increasing £. Using the distance between 
the local minimum of KE(£,r) and r=1.0 as the turbulent boundary- 
layer thickness, we see that the thickness is 

fairly constant for £>6. The turbulent boundary- layer is due 

to our eddy viscosity. The viscous sublayer, with a thickness of 
- 1/2 

KE(£,r = l) , is unresolvable . For i = 12 , the local minimum in 

KS( 5 ,, r) denotes that the turbulent boundary- layer has disappeared. 

The disappearance is because our eddy viscosity dominates the 

^. = 12 mode not just in the boundary-layer but for all values of r. 

The time-averaged fraction of the total convective flux at 

r = 0.S that is carried by modes of wavenumber S, , F^^^ ( S. ) /Total F^^^, 

con con 

is plotted in Figure 6. Far from being a smooth power law spec- 
trum, we see that: (1) over 90% of the convective flux is carried 
by the Z = 1 modes, (2) virtually no flux is carried by modes with 
.2,-6, and (3) the Jl = 3 modes advect energy downward . The convective 

flux carried by the modes of wavenumber Z, is 

Z <(e *v) „ o > • We therefore 

_ '^®r - Y,^,m -y,Z,m 

I j'*' 

expect F (iO /Total F to be proportional to the sauare root 
con con 


of the product of the kinetic and thermal energy spectra. We also 


expect 'to proportional to the correlation between the 

radial velocity and temperature: 




Y,m 


Y,m 




(3.1 


1/2 


Y,m 


We have found that there is a slow, systematic decrease in 6(1. ,r) 
with increasing £ . The decrease means that the large scale velo- 
city and temperature have better spatial correlation than the 
small scale velocity and temperature. The time-averaged correla- 
tion, <<6(£,r)>>j decreases more rapidly with £ then does <6(£,r)> 
indicating that the large scale velocity and 'temperature modes are 
also better correlated in time than are the small scale modes. No 

that «6(£=3, r=0.5» is less than zero because the £=3 convective flux is 
negative. The time-averaged total correlation of the radial velo- 
city and temperature for all £ is 

f\j- r\, 

<<6(r)>> = <<[<v^T>/ (<v^^xT^>)^^^]>> (3.2) 

In Figure 7 we have plotted <<6(r)>> as a function of radius. 

The correlation varies between 0.5 and 1.0. In comparison, 
Deardorff and Willis (1967) find that in plane-parallel convec- 

7 

tion in air (P^~ .7) with a Rayleigh number of 10 , the corre- 

lation 6(r) varies between 0.5 and 0.7. They also find a decrease 
in the correlation, 6, with wavelength. 
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C) Isotropy 

As in Paper II, we have used I(A,r)=KEpj(£ ,r)/2KEj^(A ,r) as a 
measure of the isotropy at the flow where KEj^(S,,r) is the kinetic 
energy in the radial direction at r due to motions of wavenumber Si , 
and where KEj^(£,r) is the kinetic energy in the horizontal direction. 

When lC£,r):;i, the flow is isotropic. In Paper II, it was shown 
that for all £, I(£,r) is constrained to be equal to 1 at r=0. 

Furthermore , in the outer boundary layer , the radial component of 
the velocity goes to zero, so I(£,r) goes to infinity. 

The ratio, KEpjC£ ,r) /2KEj^(£ ,r), is a useful measure of isotropy only 
for 0.2~r-0. 8. The isotropy is plotted as a function of r for all 
£ in Figure 8. The £=l isotropy is qualitatively similar to the 
values calculated with = 10^, P^ = 10 (Paper II); the radial component 

4 

of the velocity dominates except in the boundary-layer. For R =10 
we found that the radial velocity was dominant for all £, Figure 
8 shows that for £i3 the horizontal component of the velocity 
becomes important at some values of the radius. More importantly, 
for large values of £,Cl(£,r)-l] becomes small for r<0.8, 
indicating that the flow becomes isotropic at large wavenumbers . 

Note that I(£,r) for £=12 increases monotonically from r=.2 to 
ra.8; this behavior is due to the eddy viscosity. 

D) Energy Budgets 

We find that the time-averaged rate at which fluctuating 
thermal energy enters the fluid due to all modes, TE^^, is 0.166 
and the time-averaged rate at which kinetic energy enters the 
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10 


fluid due to all modes, is 6.01x10 . In Paper II we 


showed that KE.„ for an isothermal fluid is 1.8 927 rxp r £i5,9xio^°: 
in P s ’ 


we see that the calculated and isothermal values of KE . are nearly 

in 


the same. In Figure 9 we have plotted KE^^(£) and 


TEin(A), the rates at which kinetic and fluctuating 


thermal energy obtained directly from buoyancy (not from cas- 
cading) enter the shell of modes with wavenumber 1 . The defi- 
nitions TE^^(Jl) and KE^^(5,) in terms of modes appear in Paper II. 


Figure 9 shows that ^ 87% of KE . and ^ 74% of TE . are due to the 

in in 


£ = 1 modes. Note that KE£^(Ji=3) is negative,, which means that the 


5.= 3 velocity and temperature fluctuations are negatively corre- 


lated [i.e. 6(£=3,r)<0]. For £>3,KE. (£)/KE. and TE . (£)/TE. 

in in in in 


are both much less than 1 so the production part 
of the spectrum is confined to £$3. 

The rate at which fluctuating energy cascades into the set 
of modes with wavenumber £ is 


y .m JO 




TE(£). = -4tt 
cas Yj^i 

trarily write this rate as the difference 


TE C£) = TE (£) - TE, (£) 
cas up down 


(3.3) 


where 


: (£) = 

-4tt ) 

up 



Y,m 

(£) 

down 



Y,m 


-4tt T r <(v VT^), 


>r dr 


(3.4) 

(3.5) 


and where 
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= I 

l'>l 


and 

T< = T T C3.7) 

Z'<Z 

Physically, we can think of TE^(Z) as- the rate at which the ther- 

up 

mal energy cascades into the modes with wavenumber Z due to 
thermal fluctuations larger than the wavelength of the Jl-modes, and 
‘^down rate at which the thermal energy of the modes of 

wavenumber feeds the smaller thermal fluctuations. If there were 
no dissipation, then: 

TE. (£) + (3.8) 

in up down 

We plot in Figure 10 and use this ratio 

to determine whether or not the wavenumber Z belongs in the pro- 
duction or inertial range. For Jl = 0, TE^p(£) must equal zero. 

As Z increases the ratio goes smoothly to 1. For TE^^C £) / 

TEdo^n*'^^ ~ 1, which indicates that TE^^(5,)<<TE^^. (Z) or equiva- 
lently that the modes with. wavenumber Z derive most of their energy 
from the cascade and not from direct production; Z lies in the inertial 

range for Jlis. We can define KE (£) and KE,^ CZ) in a fashion 

° up down 

similar to equations (3.4) and (3,5). The ratio KE (il)/KE, ()U 

up down 

shown in Figure 10 is near unity for £>6 , indicating that the 
kinetic energy spectrum is also inertial for i~6. For A=3 the 
ratio of KE^p(Ji) to ^^d.own*'^^ anomalously large because KE^^(ii=3) 
is negative (see Figure 9). 






t'l 


Because we are numer-ically limited to a finite resolution, 

we can only approximate nonlinear quantities such as TE (H) and 

up 

TEd^wn^'^^’ approximation used in evaluating TE^ {%) is 




(3.9) 


where v is the numerically resolvable part of v (see §2). Using 
the eddy diffusivity, that was introduced in §2, we com- 

pute TE^^^^(£) with the approximation 

f r^dr ,m^ 

-'O 


J o 


dr (<T . (vVT. ) „ > 


Y,5,,m 'Se 'Y,^,m 


(3.10) 


where is the resolvable part of T^. The approximation in 
(3.11) is based on both the definition of and the assump- 

tion that the energy cascade is local in wavenumber space. We 

see from (3.11) that two terms contribute to TE, (£): one due 

down 

to the resolvable part of the temperature and the other to the 
unresolvable part that is parameterized by the eddy diffusivity. 

In Figure 11 we plot the fraction of '^^^own^''^^ that is due 
to the eddy diffusivity, TE^^^^(il ) /TE^^^^ ( £) , as a function of Z. 
Figure 11 shows that the eddy diffusivity has a negligible contri- 
bution to Z<9. This means that the eddy diffusivity 

has little effect on the large scale modes; it directly affects 
only the small scale modes . The large scale modes feel 'the 
eddy diffusivity only through the nonlinear interactions with the 
small scale modes. In Figure 11 we have also plotted the ratio. 


fraction of due to 

the unresolvable part of v parameterized through the eddy vis- 
cosities. Again, it is only the -small scale modes that directly 
feel the eddy viscosities. 

IV. TIME-DEPENDENT RESULTS 
A) Instantaneous Velocity and Temperature 

The velocity and temperature are time -dependent and are 
self-correlated for approximately one eddy 

turn-over time. We estimate the turnover time of the largest 

eddies as 2 x<<^ 2 j ^j,-1/2^2xxo~^ . Whtn averaged over many turn- 

RIO 

over times, V and T are zero. The Y ’ ’ components of the temp- 
erature, poloidal,and toroidal velocities at one particular instant 
in time are plotted in Figures 12, 13 and 14 respectively. The 
plotted functional forms of the Z-1 temperature and poloidal velo- 
city are characteristic of the forms for all values in time; it 
is only the amplitude of the functional form of the largest £=1 
mode that changes rapidly in time as the three £=1 modes take turns 
dominating the temperature and poloidal velocity. These functional 

forms are similar to those found in Paper II 

4 . ‘ . 

for R =10 , P =10. The poloidal velocity smoothly decreases from 
its maximum value at the origin to zero at the surface, with no 
nodes and with no trace of the r=0.3 shell that marks the bound- 
ary of the internal heat source. In contrast, the functional form 
of the high £ components of the poloidal velocity does change in 
time. The number of nodes in the poloidal velocity (which corre- 
sponds physically to the number of radially stacked convective 
eddies) increases with £ and is 4 for £ = 12. The functional 
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form of the a=l temperature shown in Figure 12 is doubly peaked 

with a small peak near the boundary-layer and a large peak at 

r=0.3. The toroidal component of the velocity is plotted as 
3 R 1 0 

r ’ ’ in Figure 14 and is the angular momentum about the 

z-axis. Conservation of angular momentum requires that 

I ^dr=0 for all time. The functional form of 

o 

changes dramatically in one eddy turnover time. "Jo rotation law 
can be ascribed to the flow because for all values of radius the 
time-averaged value of t \p ’ ’ is zero. It is important to realize 
that this numerical result, does not immediately follow from angular 
mon.antum conservation. 


B) Intermittency 

We have calculated some of the standard measures of inter- 

mxttency, such as the departure of the inertial spectrum from k 

and the kurtosis of the high-order structure functions. However, 

due to our limited spatial resolution, our numerically computed 

values are not quantitatively usef\^I measures of the spatial 

intermittency. Siggia and Patterson (19 7 9), working in a finely r^-solved 

Cartesian geometry have calculated the spatial intermittency for 

an isothermal fluid that is driven n'-t by buoyancy but by a 

large-scale, isotropic, permanently frozen- in velocity field. 

Instead of trying to reproduce their results in a sphere of con- 

vecting fluid, we present some quantitative measures of the time- 

dependence of the intermittency. Figure 15 shows the variations 

in the convective flux, [F -<<F >>]/<<F >> at r=0.5 as a 

con con con 

function of time. The dashed line corresponds to the variance, 
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cr, (calculated over 30 eddj^ turnover times) in units of <<F >>. 

con 

. . -7 

Each timestep plotted in Figure 15 is 3x10 units of dimension- 
less time. Each 1000 timestepsis roughly one eddy turnover time. 
There is a burst in the convective flux near iteration number 
2500. The fluctuation is of order 2a and lasts for approximately 
one half of an eddy turnover time. We have found 4 bursts similar 
to the one illustrated in Figure 15 .while examining the convective 
flux for 30 eddy turnover times. The bursts in flux are asso- 
ciated with bursts in the kinetic and thermal energy spectra. 

The fluctuations in these spectra at r = 0 . 5 

for odd values of £ are shown in Figures 16 and 17 as a function of 
time. The dashed lines correspond to the variances of the energy. 

All bursts are of order 2a and clearly cascade from the 
low to the high wavenumber modes . It takes roughly one half of an 
eddy turnover time for the' kinetic energy bursts to cascade from 
£=1 to £=9. The thermal energy cascade is slower and it requires 
nearly one turnover time to cascade from £=1 to £=9. Each of the 
bursts is initially positive, then followed by damped oscillations. 
We have never observed a burst that is initially negative. 

The fluctuations seem to arise from a secondary convective 
instability in the boundary-layer where the local Rayleigh num- 
ber is never far from its critical value. The fluctuations prop- 
agate inward from the surface via the £=1 modes, then cascade 
downward in wavelength. The fluctuations in the high wavenumber 
modes at r= . S are not caused by high wavenumber fluctuations that 
travel directly from the boundary -layer to r= . 5 , but are due to 
energy cascading from low to high wavenumbers at r= . 5 . To prove 
this, we have shown that for £>2, the 2,a fluctuation in TE(£,r=.5) 


z'-' 
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n: 
I *. 

I 


occurs just after a strong fluctuation in / . (v VT^)... 

Y r 3 I 5 " s-i 

The latter fluctuation is indicative of a cascade from big to small 


Y sJi 3in'^Y5^ sin 


modes. Furthermore, we repeated the calculation starting with 
iteration ji^lSOO as the initial data. In the repeated calculation 


we artificially set the (vVT ) . contribution to (vVT) „ 

Ys^s^l — Ys*3™ 


equal to zero. We then found that therf was still a 2a fluctuation 
in the £-1 modes, but that the fluctuations in KE(£,r=.5) 
for £>2 were less than one a. This proves that the 
£>2 fluctuations are caused by energy cascading from big to 


small modes via the (v*VT^) terms 


V. CONCLUSIONS 


We have been able to show that for sufficiently high Ray- 
leigh numbers, a convecting fluid develops an inertial subrange 


(in which the flow is nearly isotropic and the energy spectrum 
.-5/3 


= k and that part of this subrange can be resolved numeri- 

cally. The anisotropy of the flow and fine structure of the 
energy spectra are confined to the largest wavelengths. The 
effects of the unresolvable part of the inertial subrange can 
be modeled with an eddy viscosity and diffusivity that are func- 
tions only of the local values of the velocity field. The ampli- 
tudes of the eddy viscosities and diffusivity can be determined 
by requiring that the computed solutions be independent of the 
limit of spatial resolution. This method of computing the eddy 
viscosity fails to work if the limit of spatial resolution does 
not lie in the inertial subrange. Because we need a well-defined 
inertial range, our method of computing eddy viscosities requires 


ii 
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that the Reynolds and Peclet numbers are large. We have 
found that the poloidal and toroidal parts of the eddy viscosity 
are nearly equal and that the eddy diffusivity is 10% greater 
than the poloidal part of the eddy viscosity. 

The important results obtained with these eddy viscosities 
with a Rayleigh number of 10^*^ and a Prandtl number of unity can 
be summarized as follows : 

1) The fluid motions are time-dependent and have no spatial 
symmetries- 

2) An examination of the kinetic and fluctuating thermal 
energy spectra as a function of wavelength as well as 
a study of the manner in which energy cascades in and 
out of each wavelength shows that there is an inertial 
range for . 

3) A computation of the ratio of the horizontal to the 
radial component of the kinetic energy shows that flow 
tends toward isotropy at 2,~7. 

4) Except at the outer boundary-layer, nearly all of the heat 
flux is carried by convection. There is overshooting at 
the center and near the boundary- layer , causing the con- 
ductive flux to transport energy toward the center of the 
sphere . 

5) The variance, cr, (with respect to time) of the flux 
increases from the center to the surface of the sphere. 

At the surface, a is 18% of the total flux. The increase 
in a near the surface is probably due to a secondary con- 
vective instability in the boundary-layer. 
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6) The boundary-layer instabilities result in a series of 
intermittent bursts in flux and energy at all radii. 

The bursts in the energy cascade from the large, £=1, to 
the small, 1-9, wavelengths requires approximately one 
half of an eddy turnover time for the kinetic energy and 
about one eddy turnover time for the thermal energy. 
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FIGURE CAPTIONS 

Figure 1 - The time-averaged mean temperature • as a function of 
radius. The break in the curve at r=0.3 marks the bound- 
ary of the internal heat source. 

Figure 2 - The ratio of the time-averaged convective flux carried 

by all of the modes to the steady-state flux due to the inter- 
nal heat source as a function of radius. Overshooting occurs 
near the origin and boundary-layer. The variances of the 
ratio of fluxes is plotted with error bars and reaches a 
maximum near the boundary-layer. 

Figure 3 - The time-averaged 2-dimensional spectra of kinetic 
and fluctuating thermal energy as a function of 5, at r = 0.5. 

For £^7, the spectra behave as power laws. 

Figure 4 - The time-averaged kinetic energy spectrum as a func- 
tion of the 3-dimensional wavenumber, k. The dashed line 
has a slope of -5/3 and corresponds to a Kolmogorov spec- 
trum. 

2 

Figure 5 - The time-averaged kinetic energy spectra KE(Ji,r)/r 

as functions of radius for £=1,12. The spectra are similar 

u 

to the results obtained with R =10 , P =10 in that the radius 

s ^ 

at which the peak of each spectrum occurs increases with £. 

2 

Each spectrum xs normalized by xts maximum value, KE /r . 

^ ’ max 

Figure 6 - The ratio of the convective flux carried by modes of 
wavenumber £ to the total convective flux carried by all 
modes. The figure is computed at r=0.5. The negative value 






of the ratio at i-3 indicates that the 1=3 modes are per- 
forming mechanical work on the fluid. 

Figure 7 - The time-averaged correlation of the radial velocity 
and temperature <<6>> = << <v^T>/<v^>^'^ ^>> . The 

variance of the correlation is denoted by error bars. 

Figure 8 - The ratio of the horizontal kinetic energy to twice 
the radial kinetic energy as a function of radius for each 
wavenumber i. By definition, the ratio is 1 at the origin 
and goes to infinity in the boundary- layer where the radial 
component of the velocity vanishes. The ratio is a good 
measure of the isotropy for 0.2~r^0.8. As 1 increases, the 
ratio goes to unity, demonstrating that the flow becomes 
nearly isotropic at small wavelengths. 

Figure 9 - The fraction of the total amount of kinetic (thermal) 
energy generated that is produced by modes of wavenumber I 
is shown by solid (open) circles. Negligible energy pro- 
duction occurs when £~6. 

Figure 10 - The ratio of the rate at which kinetic (thermal) 

energy cascades into modes of wavenumber £ from all larger 
wavelength modes to the rate at which kinetic (thermal) energy 
cascades from the modes of wavenumber £ to all smaller wave- 
length modes. By definition the ratio must be zero for £=1. 

For £>6 the ratio approaches unity and is indicative of an 
inertial range where production and dissipation are unimportant. 


-31- 


Figure 11 - The rate at which kinetic (thermal) energy cascades 
from modes of wavenumber IL to wavelengths smaller than the 
limit of numerical resolution to the rate at which the energy 
cascades from H to all smaller wavelength modes. The rate 
at which the energy cascades to the modes smaller than the 
resolution limit is computed with the eddy viscosity Cdiffu- 
sivity) . The plotted ratio shows that only modes with i>9 
directly feel the eddy viscosity (diffusivity) . 

R 1 0 

Figure 12 - The instantaneous temperature of the Y ’ ’ mode 

as a function of radius. The functional form of this mode 
with one peak at r=0.3 and the other at the boundary-layer 
is typical of the largest £=1 mode. 


Figur’e 13 - The same as Figure 12, except the radial 

component of the velocity ,£( £+1) ^ ^/r, is shown 

here . 

Figure 14 - The same as Figure 13, except with the tor’oidal com- 
ponent of the velocity. The angular momentum about the 

3 ■ 4 

s-axis is 4Trr g j so / ^ '^r i o~^’ 

3 "^o 

functional form of r ^ g changes in an eddy turnover 
time . 

Figure 15 - The deviations of the convective flux from its time- 

averaged value <<F >> as a function of time. The. figure is 

computed at r=0.5. The dashed lines are the variance, a, 

of F in units of <<F >> . One eddy-turnover time is 

con con 

approximately equal to 1000 timesteps. 


Figure 16 -The same as Figure 15 for kinetic energy. The burst 
in energy requires ^ half a turnover time to cascade from 


il = l to il = 9. 

Figure 17 - The same as Figure 17 for the fluctuating thermal 
energy. The thermal energy cascade is slower than the 
kinetic energy cascade and requires nearly one full turn- 
over time to cascade from il = l to £ = 9. 
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